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We consider the two flavor version of the Linear Sigma Model as well as of the Nambu Jona-Lasinio 
model, at finite temperature and quark chemical potential, beyond the Mean Field Approximation. 
Using parameter values for the pion and quark current masses which weakly break chiral symmetry 
we show that both models can present more than one critical end point. In particular, we explicitly 
show that the appearance of a new critical point associated with a first order line at high temperature 
and low densities could help to conciliate some lattice results with model predictions. Using different 
techniques, we perform an extensive thermodynamical analysis to understand the physical nature 
of the different critical points. For both models, our results suggest that the new first order line 
which starts at vanishing chemical potential has a more chiral character than the usual line which 
£SJ , displays a character more reminiscent of a liquid-gas phase transition. 

PACS numbers: ll.10.Wx, 12.38.Aw, 12.39.Fe, 12.39.Ki 

(S\ ; I- INTRODUCTION 

Numerical analyses of Quantum Chromodynamics (QCD) on a discrete space-time lattice (lattice QCD), indicate 
that the transition from confined to deconfined matter at finite temperature T and vanishing quark chemical potential 
H is a crossover [1]. On the other hand, model studies [2-5] predict a first-order transition to occur for fj, of the order 
of 1/3 of the baryon mass and T = 0. In between these two regimes, a second-order critical point is expected in 
the T — fj, plane at some intermediate values of T and \x. The existence and the exact location of the critical point 

*— ', is still a matter of dispute [6] and has been under intense theoretical study using effective field theory models of 
QCD [2, 3, 7-13] (see also the recent analysis performed in Ref. [14]). Unfortunately, a direct application of lattice 
QCD at finite fi is, at present, still quite problematic. Only relatively recently, new theoretical developments and 
technical improvements allowed to circumvent in various ways the fermion determinant problem and start performing 

pg | Monte-Carlo calculations (see Ref. [15] for a review). Although most of the results obtained up to now seem to 
support the QCD critical point, an interesting observation against its existence comes from Refs. [16] where, from 
numerical simulations of QCD at imaginary chemical potential, one observes that the region of quark masses where 

r — . the transition is presumably of the first order (for quark masses smaller than the physical ones), tends to shrink for 
small positive values of the chemical potential, /z. Conversely, according to models supporting the critical point, the 
first order region should expand when /i increases, so that the physical quark mass point hits the critical line at some 
finite value of T and /i. A possible explanation for this discordance has been given in Ref. [11], where it was pointed 
out that a strong (repulsive) vector coupling may account for the initial shrinkage of the first order region, that would 
then start expanding again at larger values of fj,. As a result, two critical points might appear for a given range of 
(small) quark masses, as argued in Ref. [13]. However, it has been remarked [11] that this does not necessarily imply 
the existence of the QCD critical point since a too strong repulsive potential may in fact provoke the disappearance 
of the first order line (and thus of the critical point) for physical quark masses. If the vector coupling is too small, 
instead, the initial shrinkage of the first order region is not clearly seen. A recent estimate [17] of the vector coupling 
from flavor susceptibilities evaluated with lattice QCD seems to support the latter scenario. 

Based on the analysis of the Linear Sigma Model (LcrM) with two flavor quarks, it was shown in Refs. [12, 13] 
that the inclusion of thermal fluctuations of the mesonic fields leads to the appearance of two critical points for a 
small finite vacuum pion mass, mj < 50 MeV, without the need for a vector interaction. For physical values of the 
pion mass the model predicts only one critical point as one would naively expect for QCD. Also in this case, the 
initial shrinkage of the first order region at small fi was proved to be a not uncommon feature. However, as we will 
discuss, the direct transposition of these arguments to QCD should be done with special care. In the chiral limit, in 
fact, the phase diagram of the two flavor LctM is divided in two parts by a continuous first order transition, whereas 
in QCD universality arguments [18] suggest a second order phase transition of the 0(4) universality class, at least 
at /i = 0. Very recently, it was found (see Ref. [19]) that the correct treatment of the fermion vacuum fluctuations 
(that were neglected in most LerM applications) can change the order of the transition in the chiral limit from first to 
second order, depending on the coupling constants. In this case the phase diagram of the LoM would resemble the 
Nambu- Jona-Lasinio model (NJL) one with a second order line starting at /i = and high-T and terminating at a 
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tricritical point (at intermediate T and /i) where the first order line starts ending at T = 0. 

Interestingly enough, these findings of Rcf. [12] suggest the possibility of a rich structure for the QCD phase diagram 
in a situation which is similar to the one which arises in metamagnetic systems whose phase diagram may display 
two critical points in the magnetic field versus temperature plane [20]. A multicritical point structure induced by 
charge neutrality and vector interaction has been recently discussed by Zhang and Kunihiro in the context of the 2+1 
flavor NJL model [14]. It is worth pointing out that the presence of strong magnetic fields (B ~ 10 19 G), may change 
the order of the phase transition as shown by Fraga and Mizher [21] who considered the two flavor LctM , at [x = 0, 
obtaining that the usual crossover can turn into a first order phase transition in this regime. 

The aim of the present work is to explore in some details the phase diagram of two of the most important effective 
field theory models of QCD with two quark flavors represented by the LcrM and the NJL model when small finite pion 
(and quark current) masses are considered. The thermodynamics of both models has been compared, with standard 
parametrization in Ref. [8] where the MFA has been used. As we shall see, going beyond the MFA allows for the 
appearance of more than one critical point in both models for certain parameter values. For the LcrM we will use the 
same approximation of Rcfs. [12] and will closely follow the methods therein to map the phase diagram for various 
values of the pion mass, m° . As it will be shown, by varying m%, two and also three critical points (two of them are 
actually very close to each other) may appear. As a further step, we will explore the nature of these critical points 
by analyzing the susceptibilities and the correlations of the net-quark number density, the entropy density and the 
scalar density. 

The second part of the paper will be dedicated to the NJL model, in its simplest form, which will be treated in 
the so-called Optimized Perturbation Theory (OPT). The OPT method (which also goes by different names, or has 
many variants, e.g. "delta-expansion" [22], order-dependent mapping [23], etc.) is well known for allowing evaluations 
beyond MFA due to the way it modifies ordinary perturbative expansion, giving a non trivial (non-perturbative) 
coupling dependence. Examples of successful applications include the most precise analytical values of the critical 
temperature for non interacting Bose gases [24] as well as the precise location of the tricritical point and the mixed 
liquid-gas phase within the Gross- Neveu model in 2+1 dimensions [25]. The latter illustrates how this method can be 
a powerful tool beyond standard perturbation theory, since these important effects were missed by the MFA and could 
not be precisely determined by Monte Carlo simulations [26]. The OPT version adopted here is mainly indicated to 
non-gauge theories, which (at finite temperature) require the method to be extended, e.g. by adding and subtracting 
a hard thermal loop improvement that modifies the propagators and vertices in a self-consistent way, in the so-called 
hard-thermal- loop perturbation theory (HTLpt) [27]. Regarding its use within renormalizable theories the OPT has 
just been substantially improved by its combination with renormalization group properties [28]. 

This method has been recently applied to the NJL model in the evaluation of the thermodynamical potential beyond 
MFA using standard parametrization [29] . The same type of application will be considered here, with a different set 
of parameters, in order to verify the possibility of multiple critical points as found in the LcrM. Our investigation, 
as already anticipated, shows that in the very strong coupling limit this situation (which would be missed by the 
MFA) arises. It is interesting to remark that the OPT brings 1/N C corrections to the MFA effective potential which 
are proportional to the scalar density, p s = {ifrip), as well as to the net-quark number density, p q — (?/> + 7/>), whose 
contribution to the pressure goes as —G/(2NfN c )(2pi — p 2 s ) where G is the usual NJL coupling, A^ c is the number of 
colors and Nf the number of flavors. This means that a type of 1/N C suppressed vector term whose strength is twice 
its scalar counterpart will contribute to the pressure so that when the interaction is sufficiently strong the results seem 
to support the findings of Rcf. [11] where the SU(3) NJL with an explicit repulsive vector interaction, such as the one 
suggested in Ref. [30], was used. We recall that, contrary to the LcrM, the NJL is a non renormalizable theory which 
is often regularized by a non covariant ultra violet cut off 1 , A. Then, from the quantitative point of view our whole 
NJL application must be taken with care since to generate exotic phase diagrams similar to the LcrM one needs very 
high values for G which in turn generate high effective quark masses at zero temperature and density. Although the 
effective quark mass value generated by those parameters becomes larger than A, the values of relevant observables, 
such as the quark condensate and the pion decay constant, remain well within reasonable values. Also, as already 
emphasized, one of the goals of the present NJL application is to check whether this model, like the LcrM, supports 
the existence of more than one critical point in its phase diagram. After obtaining the OPT effective potential (or free 
energy density) we derive the pressure and many thermodynamical quantities of interest including susceptibilities and 
critical exponents. The results obtained with both models indicate that the first order line observed at low chemical 
potential and high temperature has a more "crural" character while its low temperature and high chemical potential 
counterpart displays characteristics typical of a "liquid-gas" phase transition. Finally, it will be shown that with both 
models our results seem to support the back-bending of the critical line in the \x — m c plane (where m c is the quark 



1 See Rcf. [8] for an interesting discussion regarding how this difference may affect thermodynamical results. 



current mass) which, as first discussed in Refs. [11, 13], could reconcile the actual lattice findings with most models 
predictions. It is also shown that the MFA completely misses the possibility of more than one critical point in the 
T — jjl plane, for both models. 

The paper is organized as follows. In the next Section the LcrM is reviewed with the inclusion of thermal fluctuations. 
After obtaining the phase diagram in the T — /i plane for small pion masses the characteristics of each different critical 
point is examined by a careful analysis of the densities and susceptibilities. In Sec. Ill the recent OPT application 
[29] to the NJL model is quickly reviewed. We then obtain a set of parameter values which leads to the emergence 
of a second critical point, analogous to the one found in the LcrM. We perform a comprehensive thermodynamical 
analysis to investigate how multiple critical points appear in the T — /i plane as well as in phase coexistence diagrams 
in the T — pg and P — l//9g planes. We also numerically investigate the behavior of quantities such as the interaction 
measure, the equation of state parameter, bulk viscosity, and susceptibilities at, and around, each critical end point. 
The latter quantities allow us to estimate some relevant critical exponents in order to distinguish the physical nature 
of both critical points which is also done by considering the free energy in terms of two ordering densities, p a and p q . 
Finally, in Sec IV we present our main conclusions. 

II. THE LINEAR SIGMA MODEL WITH QUARKS 

In standard notation, the density Lagrangian of the LcrM with quarks reads 



2 v " p '" j 2 



£ = o (fa) + 77 (fa) -U{<j^) + ^ [i7 M ^ - g (a + i lb r ■ tt)] V , (2.1) 



where ip is the flavor isodoublet spinor representing the quarks (u and d), and 

U(*,*) = ±(o* + **-f*) 2 -H<r, (2.2) 

is the classical potential energy density. In the chiral limit (obtained by setting H = in the previous equation) 
the chiral symmetry SU(2)v x SU(2)a is spontaneously broken at the classical level, and the pion is the associated 
massless Goldston boson. For H =^= 0, the chiral symmetry is explicitly broken by the term Ha in Eq. (2.2) which 
gives the pion a finite mass at T = and /i = 0. The scalar field a has a finite vacuum expectation value v determined 
by the classical equation of motion 

A (v 2 - f) v-H = 0. (2.3) 

Accordingly, the a field is conveniently expressed as a sum of the condensate plus fluctuations, a = v + A. In the 
LcrM Lagrangian given by Eq. (2.1) there is no explicit mass term for the quark field, the quark mass being given 
only by the condensate, gv. The parameters of the model are given by the set of equations 
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where / w = 92.4 McV is the pion decay constant. Consistently with Refs. [12, 31] we set the vacuum a mass to 
m° a = 700 MeV and the quark mass to mP — 313 MeV (i.e. one-third of the nucleon mass). The last parameter that 
one needs to fix is the vacuum pion mass, m°, which will be varied from 10 MeV to 140 McV. 

A. Linearized mesonic action for the linear sigma model with quarks 

To calculate the equation of state of the model we closely follow the self-consistent method first proposed in Ref. [31] 
and used also in Ref. [12]. Accordingly, we write the grand canonical partition function as a functional integral in the 
Euclidean space with the imaginary time r = it 



Z = Tr exp 



-13 (H - ^n) = J VipVipVa'Dir exp I f dr j d 3 x (C + pip^) \ , (2.5) 



where H and N are the Hamiltonian and the net quark number operator, respectively. In Eq. (2.5), j3 = 1/T is the 
inverse temperature, fj, is the quark chemical potential 2 and V is the system volume. Following the same steps as in 
Refs. [12, 31], we now integrate away the quark degrees of freedom. This amounts to calculate the partition function 
of the quark sector 



Zqq = / 'Di/j'Dip exp ■ 



dr / d xipDip 

/o Jv 
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where 

D = — 7°<9 r + ry ■ V — i] — g (a + vy^T ■ it) + nj . 
The Gaussian integral in Eq. (2.6) can be solved with standard techniques [32] and yields 

Z q q = det D . 
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In Ref. [33] an expression such as Eq. (2.8) has been expanded in a series of commutators involving the derivatives of 
the mesonic fields. In our analysis we will discard these terms, implicitly assuming that the meson mode amplitudes 
vary slowly in space and time (the same approximation was introduced also in Refs. [12, 31]) . The determinant in 
Eq. (2.8) can then be evaluated as for the free case. After performing a Fourier transformation in momentum-frequency 
space and using the property 
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we finally obtain 
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where Nf = 2, N c = 3, and uj n are the Matsubara frequencies taking the values w n = (2n + l)irT because of the 



antiperiodicity condition on the fermionic functional integral ip(x,0) = —ip(x,/3). In Eq. (2.10) s = Jp- 
energy, with the quark effective mass 
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We note that Eq. (2.10) is formally identical to the standard result except for the dependence of m q on the mesonic 
fields. Performing the summation over the Matsubara frequencies in Eq. (2.10) one then obtains 
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Since we are interested in the low-energy properties of the model, we will ignore (as was done in [12, 13]), for simplicity, 
the shift in the zero point energy. The effect of such a contribution at finite temperature has been analyzed in Ref. [31] 
for m v = 138 McV. For physical values of the pion mass, this is not expected to change the qualitative behavior of 
the model. As was pointed out in the introduction, however, the effect of the fermion vacuum loop, could actually 
play an important role, changing the order of the transition in the chiral limit at /i = from first to second order [19]. 
These contributions are taken into account in the NJL model where they are responsible of the dynamical breaking 
of the chiral symmetry. 

Using the result in Eq. (2.12) we can now write an effective Lagrangian that includes only the mesonic degrees of 
freedom 



- (d^ir) 2 + - [d^crf - U eff {a, it) 
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In this work wc use the quark chemical potential fi. The baryochemical potential being fig = 3fi. 
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The Euler-Lagrange equations then read: 



dd , a+ 9Ueff(^) = 0j 

da 
d^ l+ d -^l^± = 0, i = 1,2,3 
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We now proceed linearizing the mesonic action by taking the average (. . .) of the equations of motion over the held 
fluctuations. First we decompose, as before, the a field as a = v + A where v — (a) and A is the fluctuation. Of 
course, (A™) = when n is odd, and the same is true for (tt 11 ). Therefore, since the pion fluctuations always occur as 
7r 2 , the average (dU e f f / dir^ = 0, whereas from the first one of the Eqs. (2.15) we get the condition for the condensate 
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The meson masses are then identified with the average of the second derivative of the effective potential 

' d 2 U eff (v + A,tt)\ 2 / d 2 U eff (v + A, tt) ' 
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and the effective potential is linearized as 



U eff (v + A, tt) ~ {U eff (v + A, tt)) + l -ml (A 2 - (A 2 )) + \ m \ (tt 2 - (tt 2 )) 



(2.16) 
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The two terms l/2m 2 A 2 and 1/2to 2 7t 2 on the right-hand side of the last equation are the mass terms to be added to 
the kinetic energy in the mesonic Lagrangian to give the mesonic partition function 
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and 
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Stemming from Eq. (2.19), we have two self-consistency relations between the meson masses and the corresponding 
fluctuations 



(A 2 ) = 2 



dOg. 

dm 2 



on„ 



(tt 2 ) = 2^| 
x ' ami 



Finally, we write the thermodynamic potential density as 4 
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3 Also for these terms we ignore the shift in the zero point energy, as was done in [12]. 

4 In the following, where not needed, we will avoid to write explicitly the dependence of the functionals on a and 7r 2 . 



The average over the field fluctuations is performed by using the techniques of Refs. [31, 34, 35]. Given an arbitrary 
functional O (v + A, n 2 ) we can write down the Taylor expansion around A = tv 2 = and take the average term by 
term 
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By using the relation derived in Ref. [34], we then relate the terms (A n 7r 2fc ) to powers of (A 2 ) and (t 2 ), i.e. 

(A™) = (n — 1)!! (A 2 ) and (7v 2k ) = (2k + l)!!(7r 2 /3) fc , which amounts of summing up the infinite series of daisy 
and superdaisy diagrams in the Hartree approximation. The resulting expression, turns out to be equivalent to an 
integration over a Gaussian distribution [31] 



(0(v + A,tv 2 )) = J dzP rT (z)Jdyy 2 P„(y)0(v + z,y 2 ) , 



(2.25) 



where 



P<j(z) = — cxp , 

In the following, we will need to perform derivatives of the average value of some quantities such as the thermodynamic 
potential density, Q,. After two integration by parts, using Eqs. (2.25) and Eqs. (2.26) one can obtain the following 
useful relation [31] 

8, , A 2U dv /dO(v + A,ir 2 )\ 1 d (A 2 ) /d 2 (D(v + A,tv 2 )\ ld(iv 2 ) /d 2 0(v + A,ir^ 
— {0(v + A,n)) = — ^ 1 + -^—^ — 1 + -——^ _ 

With Eq. (2.27), and the Eqs. (2.16) and (2.17) it is not difficult to see that 

ki_/8u.„\ , io(A 2 ) f/ a 2 u eff \_ n ie^)(/8^ l \_ s 



dv \ dv I 2 dv \\ dA 2 / a ) 2 dv \ \ dir 
In a similar way, using the Eqs. (2.22) one can also easily show that 
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Eqs. (2.28) and (2.29) guarantee the consistency of the approach and the standard connection between thermodynamic 
and statistical mechanics [36]. The equation of state of the system is given by the simultaneous solutions in the 
variables v, m a , m^ of Eqs. (2.16) and (2.17), with the field fluctuations given by Eq. (2.29). In Ref. [12] this was done 
numerically and the authors found a rich phase structure with one or two critical points, depending on the value of 
vacuum pion mass m°. In the next subsection, we will repeat the same calculation as in Ref. [12] and we will show 
that another critical point (very difficult to detect) appears in the phase diagram. 

B. The phase diagram 

We now proceed and map the phase diagram of the model for various values of the vacuum pion mass ranging 
from m-a — 10 MeV to m v — 140 MeV. As a first step, we need to find a way to localize the first order line(s) and 
the critical end point (s). This is usually done by looking at the change in shape of the thermodynamic potential 
across a transition line. This method (that will be adopted for the NJL model in sec. Ill), allows to see quite clearly 
the onset (and usually the order) of a phase transition, but unfortunately cannot be used for the present analysis. 
Our approximation is, in fact, based on an expansion of the thermodynamic potential around a minimum and our 
equations arc defined only there so that a different method must be employed. A typical signature of the onset of a 
first order transition is the presence of two minima of the thermodynamic potential, corresponding to the high and low 
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FIG. 1: (Color online) The condensate v (normalized to /„-) as a function of the temperature in the vicinity of the first order 
transition for m n — 20 MeV and fj, = 0. The stable solution v (shown with black markers and the dashed line) is the one that 
corresponds to the highest pressure. The other two (metastable) solutions vi and Vi are also shown in the picture. 



temperature phases. In turn, we should see two distinct solutions for v,m a ,mT r at the same T and /i. The transition 
happens when these two solutions have the same potential (pressure), i.e. the point when the system switches from 
one minimum to the other. For each value of fj,, we solve the system of equations starting from low temperatures 
(the broken phase 5 ) and following the line of minima increasing the temperature by a small amount AT at each step, 
using the latest solution as initial point for the solving routine. Once we are sure to be in the high T symmetric 
phase, we solve the system of equations backwards (going from high to small T) until we are sure to be below the 
transition temperature. If the transition is continuous, we will find a unique solution for each value of T . Conversely, 
if the transition is of the first order, we will have a region where two solutions exist at the same T . 

In Fig. 1 this is shown for the condensate v around the transition line for m° = 20 MeV. For this value of the pion 
mass, the system undergoes a first order transition at T ~ 140 MeV and fi = 0. As one can see, the solution corre- 
sponding to the symmetry-broken phase (vi) exists at low T and shortly above the transition temperature, whereas 
the symmetric solution (112) exists at high T and below the transition. The actual solution (the one corresponding to 
the bigger pressure) is shown by black dots and the dashed line. 

With this method, we now map the phase diagram of the model in the T — ft plane. Our results are consistent with 
those in Ref. [12], where the authors performed the same calculation and found that for sufficiently small values of 
the vacuum pion mass (mj < 50 MeV), the system has two critical points. In Fig. 2 we plot the phase digram of the 
model in the MFA (Fig. 2a) and with the inclusion of mesonic fluctuations (Fig. 2b). In Fig. 2b we have analyzed 
various values of m° = 10, 20, 35, 50, 140 MeV. For to° = 10 MeV the phase diagram is divided in two parts by a 
continuous first order line. For mj = 20 MeV and 35 MeV the first order line is interrupted by a continuous region, 
and for some 35 < mj < 50 MeV the branch at low /j, disappears and we have the usual continuous transition at low 
fi and first order at high /i. 

The leftmost branch of the first order line (when it exists) ends in a critical point. For the rightmost one, instead, 
the situation seems to be different. A closer look revealed the presence of a double critical point. This is shown in 
Fig. 3a for mj = 35 MeV. As one can see, the first order line bifurcates and ends in two critical points that have been 
labeled C2 and c 2 . Their existence is revealed by the presence of two first order transitions that have been detected, 
as before, looking at the double solutions of our system of equations. In Fig. 3b we show the parameter v along the 
dashed line crossing the two first order lines in Fig. 3a. As one can see, there are two distinct solutions (that we have 
checked to correspond to minima of the thermodynamic potential) for two slightly different (~ 1 MeV) values of the 



5 Since the symmetry is explicitly broken by the vacuum pion mass, expressions such as "symmetric phase" or "broken phase" must be 
understood just as a nomenclature convention. 
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FIG. 2: (Color online) Left panel: the phase diagram of the LcrM in the MFA for m° = 10, 100, 140 MeV. Right panel: the 
phase diagram of the LcrM with the inclusion of mesonic fluctuations for m% = 10, 20, 35, 50, 140 MeV 



temperature, indicating a double first order. A similar behavior can be observed for all the explored values of the 
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FIG. 3: Panel a: a detail of the phase diagram for mj = 35 MeV showing the rightmost first order line bifurcating and ending 
in two critical points. Panel b: the condensate v (normalized to /„■) as a function of the temperature for fi — 248 MeV. 
The plotted region corresponds to the dashed vertical line crossing the two first order lines in Panel a. The stable solution is 
represented by black dots and the metastable solutions by open circles. 



pion mass, i.e. 20,35,50,140 MeV. Also for m^ — 10 MeV, where the transition is of the first order everywhere, a 
very short appendix of the the first order line appears around [i ~ 250 MeV and T ~ 86 MeV producing another very 
short branch ending in a critical point at \x ~ 249 MeV and T ~ 87 MeV. The situation is qualitatively similar to 
the one in Fig. 3a except for the fact that the leftmost critical point C2 does not exist as the first order line continues 
until n = 0. In Fig. 4 one can see the multiple solutions for v, indicating the two first order lines at T ~ 85.9 MeV 
and T ~ 86.2 MeV, respectively. For the rightmost transition, the discontinuity in v is very small, but its presence 
is proved by the existence of the metastable solution. At present the reason for the this unusual bifurcation of the 
high-/i first order line is not understood and requires further studies. It may simply be an artifact of the approximation 
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FIG. 4: The condensate v (normalized to /») as a function of the temperature for p = 249.6 MeV and mi — 10 MeV. Two 
first order transitions are visible at T ~ 85.9 MeV and T ~ 86.2 MeV. The stable solution is represented by black dots and the 
metastable solutions by open circles. 



(mean-field plus fluctuations) used for the treatment of the linear sigma model. However, as we will see in the next 
section, despite their vicinity in the phase digram, the two first order lines appear to have slightly different qualitative 
features. 



C. Densities and susceptibilities 



We now go one step further and analyze the susceptibilities in the LctM with the aim of characterizing the critical 
point (s) by studying the fluctuations of the net-quark number density p q , scalar density p s , and entropy density s in 
the vicinity of the critical regions. By definition, the quark density is 



T dlnZ 



Pi 



V dp 

By using the Eq. (2.27) and Eqs. (2.28) and (2.29) one then has 



on 

dp 



A, 



dU, 



eff 



dp 



T /d\nZ, 



V 



<i<i 



dp 



(2.30) 



(2.31) 



To calculate the scalar density, it is convenient to introduce a mass term —rjipip for the quarks in the Lagrangian in 
Eq. (2.1), where 77 is a fictitious bare quark mass to be set to zero afterwards. This amounts to a thermal quark mass 



f>K 



g tv + (go + nY 



From the partition function written as in Eq. (2.5) one then sees that 



Ps = (H->) = - 



TdlnZ 



V dg 



7J = 



dVL 
dr) 



r/=0 



that is, as before, 



Ps 



dg 



T /dlnZ, 



r/=0 



V 



<1<I 



dq 



(2.32) 



(2.33) 



(2.34) 



7J = 



10 

Unlike the quark density in Eq. (2.31), the scalar density can be directly evaluated, once v and the meson masses are 
known, without the need to solve any further integral. Indeed, from Eq. (2.16) one gets 



T /d\nZ a 



V 



dv 



= X(v 2 + 3(A 2 ) + (7v 2 )-f 2 )v-H 



Noting that 



dlnZ, 



q<i 



J qq ' 



drj dv 

and using Eqs. (2.34) and (2.35) we finally obtain 



di] \ dv 



1 3 In. Z 



qq 
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dv 



Ps = - g [H-X(v 2 +3(A 2 ) + ^ 2 )-f 2 )v] 



(2.35) 



(2.36) 



(2.37) 



The net-quark number density and the scalar density are the thermodynamic variables conjugate to /i and the fictitious 
quark mass rj, respectively. The last density that we need to evaluate is the the entropy density (conjugate to the 
temperature T), i.e. 



on 

"df 



In Z„ 



T dlnZ q 
V dT 



dfla dCl-n 



dT dT 



(2.38) 



In what follows, we will always assume the same number of u and d quarks in the system. The isospin density is 
therefore always vanishing for any value of T, fj, and 77, and, of course, the same is true for its derivatives with respect 
to these quantities. The Isospin density, therefore, do not mix with the other degrees of freedom 6 (p q , p s and s) and 
will be not considered in our analysis at this time. 
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FIG. 5: (Color online) Panels a and b: diagonal matrix elements of the covariance matrix as a function of the temperature 
in the vicinity of the critical point c\ for m% = 35 MeV and fi — 70 MeV and fj, — 86 MeV, respectively. Panels c and d: 
off-diagonal matrix elements of the covariance matrix as a function of the temperature in the vicinity of the critical point ci 
for mj = 35 MeV and fi = 70 MeV and [i = 86 MeV, respectively. Panels e and f: eigenvalues of the covariance matrix as 
a function of the temperature in the vicinity of the critical point c\ for 771° = 35 MeV and fi = 70 MeV and /i — 86 MeV, 
respectively. Panels g and h: square components of the (normalized) eigenvector corresponding to the larger eigenvalue (ei) of 
the covariance matrix as a function of the temperature in the vicinity of the critical point ci for m% — 35 MeV and /x = 70 MeV 
and n = 86 MeV, respectively. 



Note that even though the correlations with the other densities vanish, the Isospin susceptibility (i.e. the derivative of the Isospin with 
respect to the Isospin chemical potential) do not. 
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FIG. 6: (Color online) Panels a, b, c and d: diagonal matrix elements of the covariance matrix in the vicinity of the critical 
points C2 and c' 2 as a function of the temperature for m% — 35 MeV and p — 243 MeV,p — 246 MeV, fi = 248 MeV and 
(j, — 250 MeV, respectively. Panels e, f, g and h: off-diagonal matrix elements of the covariance matrix in the vicinity of the 
critical points ci and c' 2 as a function of the temperature for m% — 35 MeV and p — 243 MeV,^i = 246 MeV, p — 248 MeV 
and p — 250 MeV, respectively. 



We can now go ahead and calculate the 3x3 covariance matrix as 



3pq 


^Pq 


du 


dv 


Ps 


V Ps 


du 


dn 


<)s 


ds 


dji 


drj 



/-I rp I CJps Ops Ops 



,.7 = 1,2,3 



(2.39) 



This matrix is symmetric and can be diagonalized. We call ei,e2,e3 the three eigenvalues of C and u,v, z the 
corresponding eigenvectors. Each one of the three eigenvectors can be related to a different orthogonal combination 
of the three original densities p q , p s and s 



Pu = "lPq + U 2 p s + U 3 S 
Pv = VlPq+V 2 Ps +V3S 
Pz = ZlPq + Z 2 Ps + Z 3 S . 



(2.40) 



These three new densities are now independent. At the critical point, only one of them will have divergent fluctuations, 
whereas the other two will remain finite. We chose, as a convention, to order the eigenvalues from the biggest e\ to 
the smallest e 3 . At the critical point, we then expect the first eigenvalue, e 1: to diverge. According to Eq. (2.40), 
the three components of the corresponding eigenvector u = (ui, u 2 , 113) will then give us the expression of the critical 
density in terms of the original densities p q , p s and s. 

It's worth to mention that in our analysis, the correlations between p q and s are positive, whereas the correlations 
between p s and s, and p s and p q are negative. This is due to the fact that, at the transition, the thermal contribution 
to the scalar density drops from its maximum value (in the broken phase) to a very small value (exactly zero in the 
chiral limit) in the symmetric phase (see ref. [31]), whereas both s and p q exhibit the opposite behavior going from 
smaller values (in the broken phase) to higher values (in the symmetric phase). Unlike s and p q the scalar density 
is, in fact, dominated by the zero-point contributions. With their inclusion one recover the physically expected 
behavior [31]. One must then bear in mind that the following results are relevant for the thermal part of the model. 
Our general conclusions, however, are in agreement with the results obtained with the NJL in sec. Ill where the 
zero-point contributions are included. 

We analyze two different values for the vacuum pion mass: mj = 35 MeV (figs. 5,6 and 7) and mj — 140 MeV 
(figs. 8 and 9). For clarity, it is convenient to chose a label for the various critical points. We will call c\ the critical 
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FIG. 7: (Color online) Panels a,b,c and d: eigenvalues of the covariance matrix as a function of the temperature in the vicinity 
of the critical points c 2 and c 2 for m° = 35 MeV and fj, = 243 MeV,/i = 246 MeV, (X = 248 MeV and fj, = 250 MeV, 
respectively. Panels e,f,g and h: square components of the (normalized) eigenvector corresponding to the larger eigenvalue (ei) 
of the covariance matrix as a function of the temperature in the vicinity of the critical point C2 and c 2 for mj = 35 MeV and 
H = 243 MeV,/x = 246 MeV, (j, = 248 MeV and /x = 250 MeV, respectively. 
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FIG. 8: (Color online) Panels a, b, c and d: diagonal matrix elements of the covariance matrix in the vicinity of the critical 
points c 2 and c 2 as a function of the temperature for m n = 140 MeV and fi = 275 MeV,/i = 277 MeV, fi = 278 MeV and 
(j, — 280 MeV, respectively. Panels e, f, g and h: off-diagonal matrix elements of the covariance matrix in the vicinity of the 
critical points c 2 and c 2 as a function of the temperature for mj = 140 MeV and fj, = 275 MeV,/x = 277 MeV, fj, = 278 MeV 
and n = 280 MeV, respectively. 
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endpoint of the leftmost first order line (if there is any) that starts at p = and T ~ 140 MeV in the phase diagram. 
With c 2 and c 2 we will refer to the two critical endpoints of the first order line (that ultimately splits in two) that 
starts at T = and p, ~ 300 MeV, c' 2 being the rightmost of the two as in Fig. 3. As we have seen in the last 
section, for m^ = 10 MeV there is only the critical point c' 2 . We will not consider that case here, however. We 
begin with c\ for mj = 35 MeV. In Fig. 5, panels (a-b) and (c-d) we show the three diagonal components of the 
covariance matrix (Cn, C22, C33), and the three off-diagonal components (C12, C13, C23) for p = 70 MeV and 
fi = 86 MeV. For p = 70 MeV, the system undergoes a first order phase transition at T — 138.6 MeV, whereas for 
|U = 86 MeV the transition is continuous (but still very sharp) and takes place at T = 137 MeV. In between, one finds 
the unusual critical point c\ . By looking at the elements of the covariance matrix for these two values of p, one sees 
that the dominant fluctuations are given by the entropy density and the scalar density. This is shown in Fig. 5a and 
5b, where the dominant diagonal terms in the vicinity of the transition are C22 = Tdp s /drj, C33 = Tds/dT . From 
the off-diagonal terms in Fig. 5c and 5d one notices that the scalar and the entropy density are the most strongly 
anticorrelated (the magnitude of C23 = Tdp s /dT is much larger than the other off-diagonal coefficients). This fact 
can be observed in a more direct way by looking at the eigenvalue e\ and the corresponding eigenvector u. In Fig. 5e 
and 5f the three eigenvalues of the covariance matrix are plotted for the same range of values of T and p. The largest 
eigenvalue e\ is the only one showing a peak (it should actually diverge at the critical point). The eigenvector u (see 
Eq. (2.40)) gives us the three components of the critical density p u = Uip q + u 2 p s + U3S. As one can see, from Fig. 5g 
and 5h, the critical density p u is a mixture of almost solely the scalar and the entropy density (the first component 
Mi does not appear in the plot as it is always very close to zero). 
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FIG. 9: (Color online) Panels a,b,c and d: eigenvalues of the covariance matrix as a function of the temperature in the vicinity 
of the critical points c 2 and c 2 for ml = 140 MeV and p = 275 MeV, p = 277 MeV, p = 278 MeV and p = 280 MeV, 
respectively. Panels e,f,g and h: square components of the (normalized) eigenvector corresponding to the larger eigenvalue (ei) 
of the covariance matrix as a function of the temperature in the vicinity of the critical point C2 and c' 2 for mj = 140 MeV and 
p = 275 MeV,u = 277 MeV, p = 278 MeV and p = 280 MeV, respectively. 



As we have already discussed, in addition to this critical point the system exhibits two more critical points (c2 
and c 2 ) in the high-/^ region (see Fig. 3). The analysis of these two points is shown in Fig. 6 and Fig. 7 for p = 
243, 246, 248, 250 MeV. For p = 243 MeV (panels a and e in Fig. 6 and 7) the transition is continuous, but one can 
already clearly distinguish that the two peaks associated with the two critical points have already developed. The first 
one on the left (the one corresponding to c%) is very sharp, while the other (a little bit on the right, corresponding to 
c' 2 ) looks still like a bump. For p — 246 MeV (panels b and f in Fig. 6 and 7) we have a first order and a continuous 
transition. The first peak on the left is now a discontinuity, whereas the bump becomes sharper as the critical point c 2 
is approached. For p = 248 MeV we have a double first order line (panels c and g in Fig. 6 and 7) and for p — 250 MeV 
the two first order lines have merged together to form a single transition line. 

Looking at the diagonal terms of the covariance matrix in the upper panels of Fig. 6 one immediately sees that the 
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density p q is now playing a relevant role and its fluctuations are of the same order of those of the scalar density. The 
off-diagonal terms (lower panels in Fig. 6) show a positive correlation between p q and s and a negative correlation 
between p s and s, and p s and p q . 

The two critical points ci and c 2 , indeed, show different features. From the analysis of the diagonal matrix, one sees 
(lower panels in Fig. 7) that the critical density p u is in both cases dominated by the entropy density. However, For 
the point C2 (see Fig. 7 panels a and e) the scalar component is rather small (Fig. 7e), and overcomes the net-quark 
number density component only at a temperature immediately higher than the one corresponding to the peak of the 
eigenvalue e\ (compare to Fig. 7a). For the critical point c 2 , instead, the critical density p u is ~ 50% entropy density 
and the remnant 50% is almost equally split between net-quark and scalar density (see Fig. 7f). 

The same analysis has been repeated for to" = 140 McV in figs. 8 and 9 in the vicinity of the two critical points 
c 2 and c 2 (the point C\ does not exist) for p, — 275, 277, 278, 280 MeV. Once again, for p = 275 MeV the transition 
is continuous, whereas for p = 277 MeV we have a first order and a continuous transition (corresponding to the 
rightmost spike), for p = 278 MeV we have a double first order line and finally, for p = 280 MeV the two first order 
lines have merged together to form a single one. By looking at the diagonal terms of the covariance matrix (upper 
panels in Fig. 8) one can see that now the fluctuations of both critical points are dominated by the net-quark density 
and the entropy density. In addition, one notices that in all the lower panels of Fig. 8, the correlation of the scalar 
density with the entropy density and the net-quark density are almost identical. Even though in this case the scalar 
density plays only a minor role, it behaves in an even more visibly different way in the two critical points c 2 and 
c' 2 . In Fig. 8 (panels a and e) at T = 75.8 MeV we are very close to the critical point c 2 . In correspondence to this 
temperature, the scalar fluctuations C 22 are suddenly suppressed. At the same time, the correlations between p s and 
Pq (C12) and the correlations between p s and s (C23) cross zero (this can be seen, with some difficulty, in Fig. 8e). 
Indeed, in correspondence to the leftmost critical point c 2 , the scalar density seems to behave as a "spectator" of the 
critical phenomenon as it does not mix with the other densities. This is not the case for the critical point c 2 . Even if 
smaller than the entropy and net-quark number density fluctuations, now the scalar density fluctuations show a peak 
and the correlations C12 and C23 are negative. 

Increasing the mass of the pion the critical points c 2 and c 2 move to the right of the phase digram towards higher 
values of the chemical potential and lower temperature. As a result, the fluctuations of the net-quark number density 
becomes more important. In contrast, the critical point ci (when it exists), has only a minor component of the 
net-quark number density, and the transition is dominated by the entropy density and the scalar density, reflecting a 
more "crural" behavior. This is in part due to the fact that the critical point c\ exists only when the vacuum pion 
mass is small, i.e. when the chiral symmetry is only slightly (explicitly) broken. 

III. THE NAMBU-JONA-LASINIO MODEL 

Let us now consider the standard version of the two flavor NJL model which is described by [37] 

C = i> [ij^ -m c ]^ + G [{^) 2 + (^ 75 r^) 2 ] , (3.1) 

where ip (a sum over flavors and color degrees of freedom is implicit) represents a flavor isodoublet (u and d type 
of quarks) A c -plct quark fields while t are isospin Pauli matrices. The Lagrangian density (3.1) is invariant under 
(global) U(2)f x SU(N C ) and, when m c — 0, the theory is also invariant under chiral U(2)l x U(2)r groups. Note 
that, as emphasized in Refs. [30, 38], the introduction of a vector interaction term of the form (-tp^ip) 2 in Eq. (3.1) 
is also allowed by the chiral symmetry and that such a term can become important at finite densities, generating a 
saturation mechanism depending on the vector coupling strength that provides better matter stability. Regarding 
analytic nonperturbative evaluations, one can consider one loop contributions dressed up by a fermionic propagator, 
whose effective mass is determined in a self-consistent way. This approximation is known under different names, 
e.g., Hartree, large- N c or mean- field approximations (MFA). To obtain the effective potential (or Landau free energy 
density) for the quarks, U e ff, it is convenient to consider the bosonized version of the NJL, which is easily obtained by 
introducing auxiliary fields (a, 7r) through a Hubbard-Stratonovich transformation. Then, to introduce the auxiliary 
bosonic fields and to render our results more suitable to compare with the large- N c approximation it is convenient to 
use G — > X/(2N C ) and to formally treat N c as a large number, which is set to the relevant value, N c — 3, at the end 
of the evaluations. One then has 

C = $ (ij^ -m c )4>- 4>(<j + ij 5 t-iv)4> - -^-(a 2 + tt 2 ). (3.2) 

The Euler-Lagrangian equations show that a — —(X/Nc^ip = —2G^ip and n = —(X/N c )'ipij^T'ip — —2Gipi , y5Ttjj. 
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A. Optimized Perturbation Theory for the NJL model 

The basic idea of the OPT method is to deform the original Lagrangian density by adding a quadratic term like 
(1 — 8)r]tpt(} to the original Lagrangian density as well as multiplying all coupling constants by 5. The new parameter 
S is just a bookkeeping label and r\ represents an arbitrary mass parameter 7 . Perturbative calculations are then 
performed in powers of the dummy parameter S which is formally treated as small and set to the original value, 5 = 1, 
at the end 8 . Therefore, the fermionic propagators is dressed by r\ which may also be viewed as an infra red regulator 
in the case of massless theories. After a physical quantity, such as U e f/, is evaluated to the order-fc and 5 set to 
the unity a residual r\ dependence remains. Then, optimal non perturbative results can be obtained by requiring 
that U^fJij) be evaluated where it is less sensitive to variations of the arbitrary mass parameter. This requirement 
translates into the criterion known as the Principle of Minimal Sensitivity (PMS) [39] 

dUeff (k) (v) 



dq 



= . (3.3) 



In general, the solution to this equation implies in self consistent relations generating a non perturbative G dependence. 
In most cases non perturbative 1/N C corrections appear already at the first non trivial order while the MFA results 
can be recovered at any time simply by considering N c — ¥ oo. Finally, note that the OPT has the same spirit as the 
Hartree and Hartree-Fock approximation in which one also adds and subtracts a mass term. However, within these 
two traditional approximations the topology of the dressing is fixed from the start: direct (tadpole) terms for Hartree 
and direct plus exchange terms for Hartree-Fock. On the other hand, within the OPT fj acquires characteristics which 
change order by order progressively incorporating direct, exchange, vertex corrections, etc. effects. To implement the 
OPT within the NJL model one follows the prescription used in Refs. [25, 29, 40] by first interpolating the original 
four-fermion version. Then, in terms of the auxiliary fields, the deformed Lagrangian density becomes 

C = i> [vy^d* - m c - S (<r + i l5 r ■ ff ) - n (1 - 8)] V> - 5-£ {a 2 + tt 2 ) , (3.4) 

which shows that the Yukawa vertices have weight d while the meson "propagators" are proportional to 1/5. One is 
then ready to perform a perturbative evaluation of Landau's free energy in powers of 8. In the a direction, the first 
non trivial order the relevant contributions are represented in Fig. 10. Then, at finite temperature and finite chemical 
potential, the order-i5 result for the NJL free energy density is (see Ref. [29] for a more detailed discussion) 



ooe 



FIG. 10: Diagrams contributing to U e ff {fj) to order 5. The thick continuous fermionic lines represent i) = m c + 1) — 8(j] — a) 
dependent terms which must be further expanded. The dashed lines represent the a propagator and the n propagator is 
represented by dashed-doted line. The first contributes with 1/N°, the second and third diagrams (of order S) contribute with 

1/JVc 



Ueff(a) = — - 2N i N c I 1 (fi,T) + 2SN f N c (rj + m c ) ( V - a) I 2 (^,T) 

+45GN f N c if (/*, T) - 2SGN { N c { V + m c flj{^ T) , (3.5) 

where we have replaced A — ► 2GN C . In the above equation we have defined, for convenience, the following basic 
relevant integrals: 



7 Note that although we have kept the same notation, the r\ used here has as completely different role than the one used previously in the 
Lo-M. 

8 Recall that within the large- N c on performs an expansion in powers of 1/N C where N c is formally treated as large but set to the original 
value (N c = 3 in our case) at the end. 
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where e 2 = p 2 + (?7 + 77i c ) 2 . Notice also that I3 only survives at /1 7^ 0. Here, we impose a sharp non-covariant cut-off, 
A, only for the vacuum term, since the finite temperature has a natural cut-off in itself specified by the temperature. 
This choice of regularization, which allows for the Stefan-Boltzmann limit to be reproduced at high temperatures, is 
sometimes preferred in the literature [11, 29, 41]. Moreover, in the present application it assures that the temperature 
integrals appearing in the LctM and in the NJL are integrated over the same momentum range. Also, as will be 
further discussed, this regularization choice appears to be a crucial condition in order for the NJL model to reproduce 
the phase diagram with two critical points. 

The divergent integrals occurring at T = and p = are 
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(3.10) 



If needed, the T — > limit of those integrals can be readily obtained (see Ref. [29]). Then, by applying the PMS 
relation to U e ff one gets 
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(3.11) 



v=v 



Notice that if one ignores the terms proportional to G (which are of order 1/N C ) the optimal result is simply 77 = a. 
In this situation, Eq. (3.5) shows that the MFA result is exactly reproduced. Since we are mainly interested in the 
thermodynamics, one basic quantity of interest is the thermodynamical potential, $7, whose relation to the the free 
energy is given by $7 = U e ff(a). The order parameter, a is determined from the gap equation generated by minimizing 
U e ff with respect to a. From Eq. (3.5) we obtain [29] 



a = 4GNiN c (r) + m c )I 2 



(3.12) 



In order to further discuss the relation of some of our results with those obtained by Fukushima [11] it is interesting to 
note, at this stage, that since a = (a) = — 2G(^) = —2Gp s one can write the optimum thermodynamical potential 
as 



n = G P 2 s 
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(3.13) 
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where SImfal (v) has the mathematical structure of a MFA like thermodynamical potential whose effective mass given 
by fj. Since fi = —P one sees that the OPT introduces a correction like —G/(N c Nf)pl which, although suppressed 
by 1/N C , is of the same form as the one considered by the MFA when a vector interaction term like —Gvi'Vj^V^tp, 
as proposed in Ref. [30], is added to the original NJL Lagrangian density. When N c — > oo, fj — —2Gp s and the MFA 
result (for the standard NJL model) with no p q dependence is recovered. 

B. Non Standard Parametrization and the Appearance of Two Critical Points 

Usually, the three NJL parameters, G, m c and A are fixed by fitting the pion decay constant, /„. = 92.4 McV, the 
quark condensate 190MeV < — {-ipip) 1 / 3 < 260 MeV as well as the pion mass, m v — 135 MeV leading, in the MFA, 
to values such as A = 664.3 MeV, GA 2 = 2.06 and m c = 5 MeV [38]. With these values one obtains satisfactory 
predictions for the quark vacuum effective mass, m q , and for the quark condensate, (ipip) given by m° = 300 MeV 
and —{-tptp) 1 / 3 = 250.8 McV. Although in general the non covariant cut-off lies within the range 500 — 700 MeV while 
m c ~ 5 MeV the coupling can be further increased producing much higher values for m? without affecting too much the 
G independent quark condensate. For example, another MFA set also presented in Ref. [38] is given by A = 568.6 MeV, 
GA 2 = 3.17 and m c = 5.1 MeV predicting mj = 600 MeV > A and (ipip} = -247.5 McV. As one can see the value of 

m q doubles while that of the observable {ipi[J) remains within the bounds 190McV < — (-tpip) 1 / 3 < 260 MeV set by sum 
rules [42] and the value — (i/j^) 1 / 3 ~ 231 MeV which corresponds to lattice estimates [43]. Regarding applications at 
finite temperature and density one sees that also the values of the critical temperature (at p = 0) and of the critical 
chemical potential (at T = 0) increase with G. In general, the size of the first order transition line, which originates at 
high ju(~ nig) and T = 0, increases with G approaching the T axis for very high coupling strengths. This observation, 
together with the LaM results that the appearance of two critical points becomes possible for m^ < 50 MeV gives us 
the hint to use the OPT at high G and small m c (since m^ oc m c ) while setting A to usual values. For consistency we 
must use the OPT two loop relations for f n and m^ recently found in Ref. [29] and which predict some deviations from 
the MFA result for the Gell-Mann-Oakes-Renner relation. Taking A = 590 MeV, and GA 2 = 3.7 with m c = 4.5 MeV 
one obtains the very reasonable values /„. ~ 92 MeV, m, ~ 135 McV, and —(ipip) 1 / 3 ~ 264 MeV. However, now G 
has an extremely high value which is reflected in the vacuum effective quark mass value , m q ~ 787 MeV. Next, 
one can keep those values of A and G decreasing m c so as to make contact with the LcrM results. For example, 
m c = 0.1, 0.28, 0.55 MeV lead to m^ = 20, 35, 50 MeV while f n and (ipip) 1 / 3 remain very stable. With this aim let us 
analyze the phase transition pattern at high T and p = in search for a first order transition taking m c = 0.1 MeV 
(with A = 590 MeV and GA 2 = 3.7). Within the OPT this choice predicts /„ ~ 92.2 MeV, m^ ~ 20 MeV and 
— (iptp) 1 / 3 ~ 263 MeV with the expected high quark mass value, m° q — 781 MeV. In principle, one could object to the 
fact that mP > A (apart from having a numerical value much higher than 1/3 of the baryonic mass). However, as 
emphasized in the introduction, the goal of our NJL investigation has a more qualitative character and, at the same 
time, it will be shown that the relevant temperature and chemical potential for our purposes of finding evidence for 
a second critical point fall well below A. Finally, note that with our parameter choice the large value of G mostly 
affects m° q while physical observables such as /„., m n , and (-iptp) 1 ' 3 remain realistic. 

With these unusual parameter values one indeed obtains a first order transition at p = and T — 0.25 GeV as shown 
by the top panel of Fig. 11 while the usual first order transition is also observed at T = and p c ~ 680 MeV < m® 
starting a line of first order transitions as shown by the bottom panel of Fig. 11 which shows the degenerate minima 
at p = A and T = 0.096 GeV. 

The differences in latent heat, Ae, between the two transitions can be further observed by looking at Fig. 12 which 
shows the normalized e/T 4 versus the dimensionless ration T/T c . The values are Ae = 0.3 x 10~ 2 GeV 4 for p = and 
T = 250 MeV and Ae = 1.72 x 10~ 2 GeV 4 for p = A and T = 96 MeV (and of course higher for T = 0, p = 680 MeV 
but we shall refrain from numerically exploring the p > A region). 

The situation can be also observed by analyzing the thermal behavior of the order parameter represented by the 
quark condensate, v = {ipip}, which is given in Fig. 13 for p = and p = A. This figure clearly indicates that the first 
order happening at high T and vanishing p = is much softer than the one happening in the reversed situation of 
small T and high p, a fact which is well illustrated by the three dimensional plot of Fig. 14 which already shows that 
these two first order transition lines are indeed separated by a cross over region at intermediate T and p. 

This fact can be more clearly appreciated by projecting the first order lines of Fig. 14 on the T — p plane as in 
Fig. 15. This phase diagram is qualitatively very similar to the one obtained in Fig. 2 for the LctM (as mentioned 
above our value m c = 0.1 leads to m, = 20 MeV) and constitutes our first important result concerning the NJL model. 
Namely, that by going beyond the MFA and choosing the parameters so as to reproduce small pion masses one may 
obtain a second critical end point in the phase diagram. In addition, the choice of regularization procedure appears to 
be crucial for the appearance of c\ in the NJL-model. If one uses a cut off also in the T dependent integrals the critical 
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FIG. 11: Normalized free energy, U eff N (a) = [U eff (a) - U eff (a)] x 1(T 7 , as a function of a for (p = 0MeV,T c = 0.25 GeV) 
(top panel) and (fi = A,T C = 0.096 GeV) (bottom panel) showing first order phase transition. 

point c\ does not emerge. However, some authors (see Ref.[41]) have recognized that a three dimensional cut off is 
only needed at zero temperature so that the presence of high momentum quarks in the T dependent Feynman loops is 
required to ensure that the entropy density scales as T 3 at high temperature. In contrast, by being renormalizable the 
LcrM always allows for high momentum quarks in the T dependent loops. Therefore, with our regularization choice 
the temperature dependent integrals within the NJL are treated in the same footing as their LcrM counterparts. A 
comprehensive discussion about how parameters and regularization affect NJL has recently been carried out by Costa 
et al. [41]. 

For our parameter values the location of c\ happens at the point (T Cl = 249.75 MeV, /i Cl = 51.16 MeV) and of C2 at 
the point (T C2 = 149.78 MeV, /z C2 = 497.55 MeV). From the quantitative point of view these values are certainly high 
as one would expect from the fact that m® is also high. Nevertheless, the phase diagram resembles the one obtained 
with LcrM and which displays more reasonable numerical values. Note also how the first order line associated with c\ 
is almost parallel to the p axis. 

The next step is to understand the physical nature of both critical first order lines. Our study of the susceptibilities 
and densities in the LcrM has revealed that the high T small p, line terminating at c\ has a more "chiral" character 
while the small T and high p line terminating at Ci has a more hydrodynamical character. One can then map the 
T — p points to form a phase coexistence diagram such as the T versus pb/ Po shown in Fig. 16 where po — 0.17fm _ 
is the normal nuclear density. The top figure shows that the coexistence region associated with the traditional c^ 
point goes from T = to T = 149.78 MeV covering very high baryonic densities. The unusual region, on the other 
hand, is very tiny going from ps = to ps — 0.85po being restricted to a very narrow temperature range as shown 
by the bottom panel. The critical points c\ and C2 are located at ps — 0.95po and ps — 4po respectively. 

At pb — the small coexistence region terminates into a point located at T — 250 MeV. Both coexistence regions 
have very different shapes, especially near the critical points, from which one may expect that the associated critical 
exponent /?, to be evaluated later, acquires different values in each situation. Let us now map the T — jjl results into 
the P — po/ Pb plane as shown by Fig. 17 which, again, clearly indicates that the region associated with c 2 in fact 
looks like the mixed liquid-gas phase appearing in P — V type of phase diagram for a van der Waals fluid (which does 
not contain an analogous of the mixed phase associated with c\). 



C. Thermodynamical Quantities Near the Critical Points 



The OPT result for Vt = —P allows us to obtain e = —P + Ts + p,p q for the NJL with the inclusion of finite 1/N C 
corrections which makes it possible to analyze the critical behavior near each one of the two critical points. With 
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FIG. 12: e/T 4 versus T/T c for fi = 0, T c = 250 MeV (top panel) and ^ = A = 590MeV,T c = 96MeV (bottom panel). The 
associated latent heat is Ae = 0.3 x 10~ 2 GeV 4 for the top panel and Ae = 1.72 x 10~ 2 GeV 4 for the bottom panel. In both 
cases GA 2 = 3.7 and m c — 0.1. 

this aim we will numerically evaluate the following thermodynamical quantities: the interaction measure (or trace 
anomaly, A), the equation of state parameter (w = -P/e), the bulk viscosity over entropy density (C/s), the quark 
number and chiral susceptibilities (% q and Xm) as well as some critical exponents. Let us start by considering the 
interaction measure 



A = 



(e~3P) 

2^4 



(3.14) 



which is plotted in Fig. 18 for regions near c\ and ci. This quantity is expected to peak near a phase transition or 
cross over and can be helpful in locating the critical line. Fig. 18 indicates that the rise near ci looks more uniform 
around C2 and happens for temperatures near T Cl which is not surprising since there is little variation in T along the 
associated first order line as already emphasized. Although the peaks associated with c^ look more pronounced one 
has to recall that A is normalized by 1/T 4 and that the C2 region is associated with lower temperatures. 

Next, let us investigate the EoS parameter, w, as represented by Fig. 19 for the two critical points. In both 
situations one observes a downward cusp at the critical temperatures (very much like in the case of the squared speed 
of sound, yj 2 = dP/de). Below the critical temperature our results show a bump which is also observed in the usual 
LcrM , Polyakov LcM and lattice studies (see Ref.[44]). For high values of T the quantity w = P/e converges to 
approximately 1/3 as expected. 

The bulk viscosity, £, is an intrinsic dynamical quantity which, however, can be expressed in terms of the static 
thermodynamical quantities derived from the free energy as [45] 



1 
9cj 







3P) 



dT T 4 



16|e | 



(3.15) 



where ujq is a scale which will be set to A, as in Ref. [29], while eo represents the vacuum part of the energy density. 
Since £ is proportional to the specific heat, C v , the bulk viscosity over entropy density, £/s, behaves as 1/Vj 2 near 
T c in this approximation and peaks at the critical end points as indeed shown by our Fig. 20 where the divergences 
at (T Cl , /i Cl ) and at (T C2 , fj, C2 ) indicate that the energy density has a sudden change at the critical points typical of 
a first order transition. For our purposes this is an interesting quantity since it has been pointed out that one can 
distinguish whether the system experiences a first order phase transition or a crossover from observables which are 
sensitive to the bulk viscosity in RHIC type of experiments. On expects that a sharp rise of bulk viscosity near a 
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FIG. 13: The dimensionless quark condensate ratio v/vo versus T/T c for fx — 0, T c = 250 MeV (top panel) and /i = A = 
590 MeV, T c = 96 MeV (bottom panel). In both cases GA 2 = 3.7 and m c = 0.1. 
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phase transition induces an instability in the hydrodynamic flow of the plasma, and this mode will blow up tearing 
the system into droplets [44] . 

Let us now examine numerically the behavior of the quark susceptibility, Xq> & s "well as of the chiral susceptibility, 
Xm, near the two critical points. These quantities are respectively given by 



Xq = 






and 



dp s 
dm. 



(3.16) 



(3.17) 



Figures 21 and 22 show \q & n d Xm respectively as functions of T for relevant values of p. As expected, one observes 
that these two quantities peak for both Ci and c 2 . However, for Xqi the magnitude of the peak associated with c 2 
seems to be much larger while for Xm the difference is not so dramatic. In principle these results could be interpreted 
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FIG. 15: Phase diagram in the T — p plane for m c = O.lMeV (m^ ~ 20MeV) and corresponds to the choice GA 2 = 3.7. 




FIG. 16: Phase coexistence diagram in the T — Pb/po plane, where ps = p/3 is the baryonic density and po = 0.17fm~ 3 is 
the nuclear matter density. Here, m c — O.lMeV predicting m w ~ 20MeV with the choice GA 2 = 3.7. The two dark regions 
denote a mixed phase with LG denoting the liquid-gas type. The bottom panel is an expanded view of phase coexistence region 
associated with the high-T first order transition line. 



as showing that the quark density, p q , plays a very minor role within c\ while p s seems to dominate that critical point. 
At the same time, the liquid-gas type of critical point, C2, seems to receive contributions from both types of density 
in accordance with our results for the LcrM . 

These findings can be further appreciated if one investigates the free energy in terms of the two ordering densities, p s 
and p q , by using the techniques of Fujii and Ohtani [46] to Legendre transform U e ff(T, p, m c ) to V e //(T, p, m c ; p s , p q ). 
The result of this type of manipulation is shown in Fig. 23 for the two critical points, Ci and C2, as well as at an 
intermediate point (T — 240 MeV, p = 150 MeV) where a cross over takes place. The contour plots displayed by Fig. 
23 indicate that c\ is indeed dominated by the scalar interaction, p s , while the quark (vector) density also plays an 
important role at c^ in accordance with the covariance matrix results for the LcrM . 

Finally, the data used in the coexistence phase diagram T — ps allow us to evaluate the critical exponent /3 defined 
as |/3q — Pq | oc \T — Te\^ where E represents c\ or C2 while pq and p~ represent the two corresponding densities 
for a given temperature. Having the quark number susceptibility allows for the numerical evaluation of the critical 
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FIG. 17: Phase diagram in the P — po/pB plane, where pB = p/3 is the baryonic density and po = 0.17 fm~ 3 is the nuclear 
matter density. Here m c — (chiral limit) and CSB represents the broken phase ( "gas" ) while CSR represents symmetric phase 
("liquid"). The long dashed line is the T — isothermal and the dark region to its left is not accessible. The short dashed 
line represents the T = 220 MeV isothermal corresponding to a cross over temperature. The large mixed phase corresponding 
to the usual first order transition line of the liquid-gas (LG) type is labeled, the dot marks C2. The thick continuous line at 
P ~ 0.45 GeV/fm 3 represents the region associated with the unusual line and the dot marks ci. 
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FIG. 18: The interaction measure, A, for the critical point c\ (top panel) and for C2 (bottom panel). Chemical potentials above 
and below the critical points values are also shown for reference. In both cases GA 2 = 3.7 and m c — 0.1. 



exponent e, Xq tx |/i — piE\~ e , for which one may define a chiral counterpart, e m , given by \ m = \m c — m\~ e " 1 . This 
procedure is illustrated by the top panel of Fig. 24 in which e is obtained by approaching the critical point in a path 
parallel to the (i axis from fi < /i Cl ( C2 )- The values e ~ 0.64 for ci is close to the one obtained in Ref. [47] while for C\ 
the value is e ~ 0.50. The bottom panel of Fig. 24 shows the same procedure for e m with the critical points being 
approached from the m "wing" of the T — \x plane, with m > m c . Interestingly enough the numerical values get 
approximately inverted, when compared to e, and one gets e m ~ 0.49 for c-i and e m ~ 0.63 for c\. A similar type 
of procedure gives (3 ~ 0.48 for c\ and (3 ~ 0.35 for c^- Having /3 and e allows us to determine 5, 7 and a using 
e = 1 — 1/(5, 7 = j3(6 — 1) and a + 2/3 + 7 = 2. The values for the remaining exponents have been obtained after 
approximating our numbers for j3 and e by the ratios shown in Table I. Note that although our numerical estimates 
are very crude, since we do not cover many orders of magnitude nor try possible different paths leading to the critical 
points, they support our previous discussion regarding the nature of the two critical points found in this work. In 
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FIG. 19: The EoS parameter w = P/e versus the temperature for the critical point c\ (continuous line, top panel) and for C2 
(continuous line, bottom panel). In both cases GA 2 = 3.7 and m c = 0.1. For ref reference we show the EoS parameter for 
M > Mc! (dashed line, top panel) and fi < /i C2 (dashed line, bottom panel). Here, GA 2 = 3.7 and m c — 0.1. 

particular, the values for (3 listed in Table I support the liquid-gas character of en [48]. For this critical point also the 
exponent e, associated with x q , is greater than the one associated with the critical point C\. Note also that the values 
of a obtained by using our approximate ratios for /3 and e in the scaling relations show that a = e for both critical 
points which is consistent with the universal arguments presented in Refs. [10, 49] since it is expected that Xq an d 
C v should be the same near the critical points. A detailed and highly precise determination of the associated critical 
exponents is beyond the scope of the present application and the interested reader is referred to Refs. [10, 41, 47]. 

TABLE I: Approximate ratios for the critical exponents associated with c\ and C2 in the NJL model. The values in between 
brackets are the results which have been numerically obtained. 
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D. Back-bending in the fi 



plane 



It will be interesting to explicitly show that, at least from a qualitative perspective, our model predictions can be 
relevant for the lattice results obtained by de Forcrand and Philipsen [16] who observed a shrinkage of the first order 
transition region when considering higher values of \i. For completeness let us also compare the OPT with the MFA 
results when the model parameters are tuned so that the latter approximation also generates a first order transition 
at p, = 0. Within the NJL model, the MFA can predict a first order transition at fi = and for A = 590 MeV when 
the higher coupling, GA 2 = 3.98, is used with m c = 0.1 MeV. In this case, the MFA also predicts f n ~ 93 MeV, 
m^ ~ 20 MeV, and -(^V>) 1/3 ^ 263.5 MeV, and m° ~ 814MeV while, for the same set of parameters, the OPT 
predicts U ^ 90.68 MeV, m^ ~ 20.73 MeV, and -(i/^) 1 / 3 ~ 265 MeV, and m° ~ 856 MeV. 

In the /i — m c plane the NJL model then generates Fig. 25 in which the top panel, corresponding to the MFA, shows 
only a one branch first order line. When the OPT is applied to the NJL with the stronger GA 2 = 3.98 considered 
in this subsection one observes a first order line for m c < 3.5 MeV which roughly corresponds to m, = 122 MeV. 
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FIG. 20: Bulk viscosity over entropy density, £/s, as a function of the temperature for the critical point c\ (top panel) and for 
C2 (bottom panel). Chemical potential above and below the critical points values are also shown for reference. Here, GA 2 = 3.7 
and m c — 0.1. 

Starting at the point m c — 3.5 MeV, ll = one may follow a second order transition line which goes left touching the 
ii axis at [i ~ 70 MeV and then going up to fi ~ 350 MeV where it bends back to the right hand side for finite values 
of m c . The situation is illustrated by the bottom panel of figure 25. 

Therefore, when one goes beyond MFA, the NJL also predicts a bending back behavior of the critical line in a 
way which is consistent with the negative curvature of the lattice simulations of Refs. [16]. This result also supports 
Fukushima's suggestion [11] regarding the eventual back bending. Although this author has considered the SU(3) 
version of the NJL (being able to reproduce a critical surface in the m u ^ — m s — \x space) with standard parametrization 
it is interesting to note that the bending was obtained by adding a vector interaction which generated a — Gypi 
contribution to the pressure. By looking at our Eq. (3.13) one sees that the OPT also brings a term like —G/{N c Nf)pl 
to P = — f2. Finally, note that a back bending behavior is also implied by our LerM results if one considers that c-2, 
and Co almost coincide. 



IV. CONCLUSIONS 



In recent publications [12, 13] it was shown that for small values of the vacuum pion mass (< 50 MeV) the phase 
diagram of the LerM has two distinct first order lines terminating in two critical points. One, as predicted by most 
models, starts at finite li and T = 0, whereas the other is an unusual first order line which starts at high T and \i = 0. 
In between the critical endpoints (whose exact location depends on the vacuum pion mass) of these two lines there is 
a crossover transition. At the origin of this behavior, that is not observed in MFA, are the thermal fluctuations of the 
mesonic fields, that have been taken into account adopting a self consistent method first proposed in [31]. Inspired 
by this finding we have performed a careful analysis of the LerM with the aim of understanding in deeper details 
the nature of the critical points. At the same time, we have considered another popular effective quark model, the 
NJL model, to investigate under which conditions a phase diagram with at least two distinct critical points could be 
reproduced. 

The analysis of the LerM has been performed by using the same method of Ref. [12] and our results agree with the 
ones obtained therein: the unusual first order line which starts at high T and ll = ending in the "new" critical point 
(ci) was also found. At the same time, a closer look in the vicinity of the "usual" critical endpoint revealed that, 
right before to end, the first order line that starts at finite ll and T — bifurcates in two critical endpoints (c2 and 
c' 2 ) separated by few MeV in T and ll. 

In addition to the phase diagram we have studied susceptibilities and correlations of the net-quark number density, 
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FIG. 21: Normalized quark susceptibility, Xq/A 2 , f° r the critical point c\ (top panel) and for C2 (bottom panel). Chemical 
potentials above and below the critical points values are also shown for ref reference. Here, GA 2 = 3.7 and m c — 0.1. 
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FIG. 22: Normalized chiral susceptibility, Xm/A 2 , for the critical point ci (top panel) and for C2 (bottom panel). Chemical 
potentials above and below the critical points values are also shown for ref reference. Here, GA 2 = 3.7 and m c — 0.1. 



the entropy density, and the scalar density. Our results for the covariant matrix show that for c\ the dominant 
fluctuations are given by the scalar density and the entropy density, whereas the quark number density plays only 
a minor role. For the two critical points C2 and c' 2 , instead, also the fluctuations of the quark number density are 
important. Their contribution becomes larger and larger as the vacuum pion mass is increased, probably also due to 
the fact that ci and c' 2 move towards higher values of [i. Despite their vicinity in the phase diagram, ci and c' 2 seem 
to exhibit distinct features. The main difference being the fluctuations of the scalar density, that are more important 
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FIG. 23: The Legendre transformed effective potential, V e //(T, p, m c ; p s , p q ), projected in the dimensionless plane p a — p q plane, 
where p = p/(10 6 MeV 3 ) in both cases. The top panel shows this quantity at c\, the middle panel shows it at T — 240 MeV, p = 
150 MeV where a cross over takes place. The bottom panel shows the contour plot at C2 
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FIG. 24: (color on line) Top panel: Logarithmic plot of the dimensionless Xq/A 2 as a function of A(p) — \fi — p C2 (c 1 ) approached 
from p < Pc 2 (c 1 )- The dots correspond to ci and the diamonds to ci. Bottom panel: Same type of plot for the dimensionless 
Xm/A 2 as a function of A(m) = \m — m c \ approached from m > m c . 
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We have then considered the NJL thermodynamical potential, recently evaluated with the OPT [29], with parameter 
values which simulate small pion masses in order to generate a T — p phase diagram where the existence of two critical 
points, c\ and C2, separated by a cross over region has been observed. This type of phase diagram is similar to the 
one found, in the temperature-magnetic field plane, for the compressible metamagnetic Ising model [20] 

We have performed an extensive thermodynamical analysis in order to find the essential physical features that 
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FIG. 25: The /j, versus m c plane showing the second order transition boundary (dashed lines) associated with the critical 
point(s) for the NJL. The shadowed regions correspond to first order transitions and the white region represents the cross over 
region. The top panel corresponds to MFT and shows only a single first order branch with one critical point. The bottom 
panel corresponds to OPT and shows the possibility of two first order branches associated with two distinct critical points. 
Both results are for GA 2 = 3.98. 



distinguish both critical points concluding, in agreement with the LcrM case, that the usual one (C2), located at 
intermediate values of T and fj, has a more hydrodynamical character with the quark number density playing an 
essential role. This density has little influence at lower values of fj, where the new, unusual, critical point c\ region 
is dominated by the scalar density. Although we have not considered a vector term in the version of the NJL model 
considered here it is rather interesting to note that the OPT pressure has the term —G/(2NfN c )(2pi — p 2 s ) which by 
being 1/N C suppressed does not contribute to the MFA which completely misses out the possible existence of c\ and 
the associated first order line. 

From these results it is clear that the inclusion of contributions beyond MFA can have dramatic consequences on 
the phase diagram of quark models and cannot be neglected. With the combined analysis of the LcrM and the NJL 
model, we have shown that the appearance of multiple critical points is not an uncommon feature of effective quark 
models, even without explicitly introducing a vector interaction as was done in [11]. Before exporting these notions 
to QCD, however, some caution is advisable. As has been very recently shown in Ref. [19], based on the analysis of 
the LcrM in MFA, the inclusion of vacuum fluctuations can change the transition in the chiral limit at /i = from first 
to second order so that this model would behave as the NJL model in the same regime. However, our results for the 
latter suggest that the consideration of vacuum contributions should not influence the appearance of the critical point 
ci, at least qualitatively, when an appropriate tuning of model parameters is carried out beyond the MFA. In the 
minimal version of the NJL model considered here ci appears only if the coupling is larger than usual (see discussion 
in sec. IIIB) while a sharp cut off is introduced only in the divergent integrals. However, having a large parameter 
space, more sophisticated versions of the model could also allow for the appearance of a critical point point like c\ 
generating, at the same time, more plausible values for the vacuum effective quark mass. Concerning the other two 
critical points found in the LcrM (ci and c 2 ), instead, vacuum fermion loops are not expected to modify dramatically 
the qualitative behavior of the model, especially for values of the vacuum pion mass close to the physical one [31]. 
The presence of c' 2 in the LcrM with thermal fluctuations, however, has not been confirmed by the OPT analysis of the 
NJL model. This might be due to various reasons including the fact that the two models have been evaluated under 
different approximations, the fact that the NJL does not have true mesonic degrees of freedom, etc. The clarification 
of these points calls for further investigations. Nevertheless, our results allow us to conclude that it is possible to use 
these effective models to generate metamagnetic like phase diagrams at the expense of using non standard parameters 
for m c and m^ which weakly break chiral symmetry in approximations which go beyond the MFA. Having two critical 
points, it then becomes possible to observe the shrinkage of first order region, in the /j, — m c (Tn n ) plane, as one 
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considers higher values of /j, so that model predictions could be conciliated with the lattice results by de Forcrand 
and Philipsen [16]. Then, as we have shown, the first order region should increase again eventually intercepting the 
physical mass point as argued in Ref.[13]. The "crural" like first order transition observed in the line associated with 
C\ seems to be much softer and the coexistence regions much smaller than in the traditional "liquid-gas" type of 
line associated with the c-i which perhaps would make it harder to be detected in lattice QCD calculations (see for 
example ref. [50]). Although quantities such as the trace anomaly, EoS parameter, bulk viscosity and susceptibilities 
display the expected behavior associated with first order phase transitions it seems that the different critical points 
belong to distinct universality classes. 
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